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ABSTRACT 


' :  1 


In  this  report,  we  demonstrate  the  use  of  adaptive  homomorphic  fil¬ 
tering  in  exposing  objects  under  light  cloud  cover.  More  specifically, 
the  homomorphic  filter  invoked  is  space-varying  and  is  parameterized  by 
the  local  mean  level  of  the  degraded  image  which  serves  as  an  indication 
of  the  extent  of  local  cloud  cover  degradation.  This  approach  repre¬ 
sents  a  departure  from  other  attempts  to  enhance  similarly  degraded  im¬ 
ages  in  that  they  have  relied  on  nonadaptive  (long-space)  homomorphic 
filtering. 

We  show  through  some  preliminary  experiments  that  the  adaptive  proce¬ 
dure  compares  favorably  with  long-space  stochastic  homomorphic  filter¬ 
ing.  In  particular,  adaptive  homomorphic  filtering  appears  to  have 
greater  potential  than  the  long-space  methods  in  exposing  objects  be¬ 
neath  light  cloud  cover,  while  it  does  equally  well  as  an  adaptive 
high-pass  filtering  technique.  In  addition,  adaptive  homomorphic  fil¬ 
tering  compares  favorably  with  an  iterative  homomorphic  enhancement 
procedure  which  we  have  devised  to  improve  results  from  the  one-pass 
nonadaptive  homomorphic  filter. 
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1 .  INTRODUCTION 


In  a  number  of  different  applied  problems  in  optical  imaging,  a  mul¬ 
tiplicative  degradation  is  often  introduced  in  the  imaging  process.  For 
example,  an  image  obscured  by  shadows  or  light  cloud  cover  can  be  mo¬ 
deled  as  the  product  of  a  function  representing  the  original  object  or 
scene  reflectivity  and  a  "noise”  component  representing  the  degradation. 

In  such  problems,  the  nonlinear  procedure  of  homomorphic  filtering 
has  been  used  successfully  in  image  restoration  and  enhancement.  In 
particular,  Oppenheim  et  al  [1]  and  Fries  and  Modestino  [2]  have  applied 
homomorphic  filtering  to  enhance  images  whose  fine  structure  has  been 
obscured  by  the  effects  of  shadows.  The  former  procedure  relies  on  a 
deterministic  model  of  shadow  formation,  while  the  latter  relies  on  a 
stochastic  model  which  Invokes  a  random  field  representation  of  shadows. 
In  both  approaches,  a  single  homomorphic  filter  is  applied  to  the  image 
to  be  enhanced;  i.e.,  one  filter  is  applied  on  a  long-space  basis. 

A  similar  long-space  (stochastic)  approach  was  taken  also  by  Mitchell 

et  al  [3]  in  restoring  images  degraded  by  light  cloud  cover.  Here  a 

particular  cloud  pattern  is  viewed  as  a  sample  function  of  a  random 

field.  Based  on  estimated  cloud  statistics,  a  fixed,  long-space  homo- 

/ 

morphic  filter  was  applied  over  the  entire  image  to  expose  objects  under 
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light  cloud  cover.  In  neither  of  these  problems  ( i . e ., restoring  images 
degraded  by  shadows  or  light  cloud  cover)  were  the  local  properties  of 
the  degraded  image  considered. 

In  this  report,  we  take  an  alternate  approach  to  the  particular  prob¬ 
lem  of  enhancing  images  degraded  by  light  cloud  cover.  As  did  Mitchell 
et  al  [3],  we  apply  a  homomorphic  filter.  This  filter,  however,  is 
adaptive  and  requires  local  deterministic  quantities;  i.e.,  the  filter 
is  space- varying  and  is  parameterized  by  the  local  mean  level  of  the  de¬ 
graded  image  which  serves  as  an  indication  of  the  extent  of  local  cloud 
cover  degradation.  This  approach  is  similar  in  style  to  the  adaptive 
image  enhancement  and  restoration  techniques  of  Peli  and  Lim  [4J,  Lira 
[5],  Schreiber  [6],  and  Gilkes  [7]. 

The  stochastic  model  of  cloudy  images  invoked  by  Mitchell  et  al  [3] 
is  effectively  the  product  of  a  function  of  the  ideal  image  and  a  cloud 
transmission  function  defined  over  the  entire  extent  of  the  image.  Our 
cloudy  image  model,  on  the  other  hand,  is  deterministic  in  nature  and 
relies  on  dividing  the  image  into  a  number  of  overlapping  windowed  seg¬ 
ments  to  each  of  which  we  associate  a  particular  short-space  cloud  tran¬ 
smission  function.  The  logarithm  of  each  windowed  segment  can  be  divid¬ 
ed  into  two  approximately  disjoint  spectral  bands:  the  cloud 
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transmission  function  being  low-pass  and  a  function  of  the  image  being 
high-pass.  Although  the  image  contains  low-pass  information,  we  shall 
assume  it  is  not  very  significant  in  exposing  object  shapes  under  light 
cloud  cover.  In  addition,  we  depart  from  a  typical  assmption  that  the 
desired  image  itself  can  be  modeled  as  the  product  of  illumination,  a 
low  frequency  component,  and  reflectivity,  a  high  frequency  component. 
Rather,  we  view  the  illumination  component  in  the  same  way  as  we  view 
the  reflectance  component,  i.e.,  as  having  an  important  high  frequency 
component  due  to  the  interaction  of  light  and  ground  objects  [  6 ] . 
Furthermore,  we  assume  this  high  frequency  component  is  approximately 
disjoint  from  the  cloud  transfer  function  which  is  low-pass. 

To  stage  a  fair  comparison  between  our  approach  and  that  of  Mitchell 
et  al  [7],  we  have  attempted  to  exhibit  the  best  of  both  techniques. 
Toward  this  end,  we  were  led  to  an  improvement  of  the  latter  algorithm. 
Specifically,  this  improvement  invokes  various  iterative  techniques  for 
obtaining  a  better  signal  estimate  by  repeatedly  updating  the  Wiener 
filter  required  in  their  algorithm.  We  have  also  included  In  our  com¬ 
parison  the  result  of  adaptive  high-pass  filtering  and  long-space  deter¬ 
ministic  homomorphic  filtering. 
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Before  proceeding  with  the  development  of  our  new  methods  and  compar¬ 
isons,  we  review  some  important  ideas  and  formulate  a  framework  for  our 


investigations. 

t 


2.  MODELING  THE  IMAGING  PROCESS 

A  number  of  approaches  to  modeling  the  imaging  process  have  been  pre¬ 
sented  in  the  literature.  A  more  difficult  problem  is  to  model  the  im¬ 
aging  process  in  the  presence  of  light  cloud  cover.  In  this  section,  we 
first  review  one  traditional  approach  to  modeling  undegraded  images. 
Two  different  viewpoints  are  presented  which  rely  on  an  illumination-re¬ 
flectivity  model.  We  then  enter  a  stochastic  framework  in  which  a  model 
of  Images  degraded  by  light  cloud  cover  will  be  discussed.  Finally,  we 
present  our  own  deterministic  interpretation  of  cloudy  images.  As  we 
shall  see  in  the  following  sections,  each  model  of  a  light  cloud  covered 
image  leads  to  a  specific  image  enhancement  procedure. 
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2.1  Deterministic  Modeling  of  Undegraded  Images 

In  the  formation  of  Images,  the  illumination  and  reflectivity  are 
combined  by  a  multiplication  law  [1]: 

I(n,m)  -  i(n,m)  r(n,m)  (1) 

where  i(n,m)  and  r(n,m)  are  the  illuminance  and  reflectance  components 
respectively.  One  assumption  is  that  the  illumination  varies  slowly 
(i.e.,  it  contains  mainly  low  frequency  components)  and  the  reflectance 
is  sometimes  dynamic  and  sometimes  static  and  therefore  may  be  regarded 
as  containing  mainly  high  frequency  components.  The  logarithm  operation 
separates  the  multiplicative  signal  into  two  additive  components: 

log[l(n,m)]  -  log[i(n,m)]  +  log[r(n,m)]  (2) 

If  a  linear  amplifier  or  attenuator  with  gain  a  follows  the  log,  the 
image  output  obtained  by  exponentiating  alog[I(n,m)}is  given  by 

I'(n,m)  -  [I(n,m)]a  (3) 

When  a<l,  we  obtain  a  washed  out  image  and  when  a  >1  ,  the  image  is 
sharpened,  but  might  be  saturated  [1].  Therefore,  if  we  want  to  reduce 
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the  dynamic  range,  which  Is  contributed  mostly  by  the  illumination, 


while  simultaneously  increase  the  sharpness,  which  is  contributed  mostly 
by  the  reflectance,  we  should  choose  a  as  a  function  of  frequency, 
i.e.,  a  <1  for  low  frequencies  and  a  >1  for  high  frequencies.  More  spe¬ 
cifically,  we  would  apply  a  two-dimensional  ( 2  — D )  high-pass  filter 
H(u,v)  to  the  2-D  Fourier  transform  of  log[I(n,m)].  The  cross  section 
of  the  shape  of  the  required  filter  is  illustrated  in  Figure  1.  Expo¬ 
nentiating  the  modified  Fourier  transform  yields  an  enhanced  image. 
This  procedure  for  image  enhancement,  termed  homomorphic  filtering,  was 
first  proposed  by  Oppenheim  et  al  [1].  Since  the  illumination  is  not 
truly  low-frequency  and  the  reflectance  consists  not  only  of  high  fre¬ 
quency  components  but  also  of  low  frequency  components,  the  output  will 
have  certain  artifacts.  It  seems  that  the  degree  to  which  those  arti¬ 
facts  are  visible  is  strongly  dependent  on  the  way  a.  changes  from  a  <1 
to  a >1  [1], 


Fig.  1.  Cross  section  of  a  filter  with  variable  a  . 


In  an  alternate  viewpoint,  Schrieber  [6]  hypothesizes  that  in  natural 
scenes  the  illumination  function  is  as  essential  to  perception  as  the 
reflectivity  of  the  object.  Both  have  a  broad  spectrum,  generally  with 
more  power  at  the  low  spatial  frequencies.  The  reflectance  derives  its 
low  frequency  component  from  the  presence  of  large  patches  of  relatively 
constant  value.  The  illumination  derives  its  high  frequency  component 
from  interaction  between  the  incident  light  and  the  edges  and  surfaces 
of  objects,  which  are  at  many  different  angles.  Consequently,  any  at¬ 
tempts  to  separate  these  components  by  long-space  homomorphic  filtering 
will  be  limited  in  their  success.  Rather,  an  adaptive  control  of  the 
low  and  high  frequency  components  at  each  point  (along  with  neighboring 
points)  of  the  original  scene  may  be  preferable  for  enhancement.  Such 
processing,  sometimes  referred  to  as  short-space  processing  is  imple¬ 
mented  in  the  system  of  Figure  2,  as  proposed  by  Schrieber  [6J.  In  this 
general  processor,  we  can  enhance  the  high  frequencies  and  reduce  the 
dynamic  range  on  a  short-space  basis  in  regions  where  it  is  needed.  An 
adaptive  process  in  this  same  spirit  could  of  course  be  applied  to  the 
logarithm  of  the  image,  resulting  in  an  adaptive  homomorphic  procedure. 
(It  is  interesting  to  note  that  Schrieber  has  shown  that  under  a  low- 
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contrast  condition,  a  single  (long-space)  homomorphic  filter  approxi¬ 
mates  a  single  linear  filter  followed  by  a  nonlinearity.)  The  sectioned 
method  that  we  propose  in  Section  3  of  this  report  represents  a  cross 
between  the  long-space  homomorphic  (with  a  single  filter)  and  Schrie- 
ber"s  short-space  adaptive  approach  (which  effectively  invokes  a  filter 
that  is  changing  at  each  pixel). 


Fig.  2.  Adaptive  control  for  image  enhancement. 


2.2  A  Stochastic  Approach  to  Modeling  Images  Degraded  by  Light  Cloud 
Cover 

In  the  model  of  Mitchell  et  al  [3],  it  is  assumed  that  some  energy 
of  the  ground  reflection  passes  through  the  clouds.  The  energy  collect¬ 
ed  above  the  clouds  consists  of  two  components:  ground  information  and 
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scattering  from  the  clouds.  Figure  3  is  a  simplified  model  of  this  im¬ 
aging  process.  The  energy  recorded  is  given  by. 

s(n,m)  »  aLr(n,m)t(n,m)  +  L[l-t(n,m)]  (4) 


where 

L  ”  sum  illumination 

t(n,m)  “  transmittance  function  of  the  cloud  (0<t(n,m)j<l )  in  a 
downward  direction 


a  ■  attenuation  of  the  illumination  in  a  downward  direction 


from  the  clouds  to  the  ground  (assumed  constant) 


r(n,m)  **  reflectance  of  the  ground 


GROUNO 

REFLECTANCE 


Fig.  3.  A  simplified  model  of  the  imaging  process  above  clouds. 
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In  order  to  recover  the  desired  reflectance  information,  it  was  as¬ 
sumed  that  the  transmittance  function  of  the  cloud  and  the  reflectance 
of  the  ground  are  stochastic  processes.  It  was  also  assumed  that  the 
transfer  function  of  the  cloud  has  relatively  more  energy  in  low  spatial 
frequencies.  From  (4),  we  can  obtain  a  multiplicative  form  of  the 
’noise"  and  the  "signal"  given  by  (5): 

L-s(n,m)  =  t(n,m)  [L-aLr(n,m)]  (5) 

By  taking  the  log  of  (5),  we  obtain  the  addition  of  the  "signal”  and 
"noise”  which  can  be  written  as 

log[L-s(n,m) ]  =  log[t(n,m)]  +  log[L-aLr(n,m) ]  (6a) 

or 

P(n,m)  =  N(n,m)  +  M(n,m)  (6b) 


where 


P(n,m)  =*  log  [L-s(n,m)]  =  "signal  +  noise" 
N(n,m)  =  log[t(n,m)]  *  "noise” 

M(n,m)  -  log[L-aLr(n,m) ]  *  "signal" 


and  where  P(n,m),  U(n,m)  and  M(n,m)  are  viewed  as  stochastic  processes. 


One  approach  to  estimating  the  signal  from  (6)  involves  the  applica¬ 
tion  of  Wiener  filtering.  If  we  assume  that  the  signal  and  the  noise 

are  uncorrelated  (a  reasonable  assumption  here  since  the  formation  of 
clouds  has  effectively  no  correlation  with  the  formation  of  the  ground 
images),  we  can  obtain  the  minimum  mean  square  error  estimate  of  the 
signal  by  filtering  the  given  degraded  image  with  the  following  optimal 
linear  filter: 

H(u,v)  =  SMp(u,v)  /  Spp(u,v) 

[Spp(u,v)  -  S^Cu.v)  -  S(u,v)]  /  Spp(u, v)  (7) 

=  the  power  spectrum  of  the  signal  plus  noise 

•  the  power  spectrum  of  the  noise 

*  the  cross  power  spectrum  between  the  signal  and  noise 

“  the  mean  value  of  the  signal 

■  the  mean  value  of  the  noise 

■  the  2-D  dirac  delta  function 


where , 


spp<u,v) 


SNN(U’V) 


sMp(u,v) 


6(u,v) 


H 


Mitchell  et  al  [3]  chose  to  estimate  Spp(u,v)  as  the  magnitude 
squared  of  the  2-D  discrete  Fourier  transform  of  log [L-S (n,m) ]  where  L 
is  estimated  as  the  highest  value  in  the  image.  S^Cu.v)  was  estimated 
through  an  estimate  of  t(n,m)  given  by 

L-s(n,m) 

t(n,m)  *  — — — —  (8) 

L-G 

where  G  is  a  constant  which  ideally  should  equal  the  typical  ground  re¬ 
flection  (i.e.,  the  average  of  aLr(n,m)).  The  power  spectrum  of  the 
noise  was  then  roughly  estimated  by  the  magnitude  squared  of  the  Fourier 
transform  of  t(n,m).  Since  the  transfer  function  of  the  cloud  transmit¬ 
tance  was  considered  to  possess  only  low  frequencies,  the  estimate  of 
this  power  spectrum  was  low-pass  filtered.  Finally,  an  estimate  of 
and  can  be  formed  from  the  initial  values  of  the  correlation  func¬ 
tions  associated  with  the  estimates  of  Spp(u,v)  and  S^Cu.v)  and  the  re¬ 
lation  given  by  (6).  The  cross  section  of  the  general  shape  of  the  re¬ 
sulting  optimal  filter  is  given  in  Figure  4. 
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Fig.  4.  Cross  section  of  the  general  shape  of  the  optimal  Wiener  fil¬ 
ter  . 

The  filter  has  a  boost  at  (u,v)=*(0,0)  as  a  result  of  assuming  a  non¬ 
zero  mean  process  and  an  attenuation  at  the  low  frequencies.  Excluding 
the  value  at  (0,0),  it  has  the  same  shape  as  that  used  in  the  homomorph¬ 
ic  filtering  procedure  of  Oppenheim  et  al  [1],  described  in  Section  2.1 
i.e.,  the  resulting  filter  for  the  stochatic  approach  leads  to  a  special 
case  of  the  filter  for  the  deterministic  homomorphic  approach.  It 
should  be  noticed  that  in  the  stochastic  approach  the  resulting  filter 
is  fixed  for  a  specific  image  and  cloud  cover,  while  in  homomorphic  fil¬ 
tering  the  values  of  the  filter  at  the  high  frequencies  are  variable 
i.e.,  we  can  increase  them  more  than  in  the  stochastic  approach. 
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2.3  A  Deterministic  Approach  to  Modeling  Cloudy  Images 

Our  approach  to  the  problem  of  modeling  cloudy  images  invokes  the  as¬ 
sumption  that  images  and  clouds  are  deterministic  processes.  First  we 
adopt  Schreiber's  approach  [6]  which  assumes  that  both  illumination  and 
reflectance  contain  low  and  high  components  and  that  the  high  frequen¬ 
cies  are  important  for  perception.  At  the  same  time,  we  assume  that  the 
transfer  function  of  the  cloud  transmittance  contains  only  low  frequen¬ 
cies  on  a  short-space  basis.  We  now  interpret  (6)  so  that  t(n,m)  (de¬ 
fined  on  a  short-space  basis)  is  deterministic  and  contains  mainly  low 
frequencies  and  L-aLr(n.m)  is  also  deterministic  and  contains  perceptu¬ 
ally  important  high  frequency  information.  By  using  these  assumptions, 
the  image  aLr(n,m)  can  be  estimated  in  a  deterministic  fashion  through 
homomorphic  filtering  which  reduces  the  low-frequency  contribution  of 
t(n,m).  Furthermore,  the  estimate  can  be  enhanced  through  homomorphic 
filtering  by  amplifying  high  frequency  components.  We  elaborate  on  this 
procedure  in  the  next  section. 
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3.  EXTENSIONS  OF  HOMOMORPHIC  FILTERING 


The  assumption  that  clouds  are  stationary  and  approximately  low-pass 
motivated  the  long-space  stochastic  homomorphic  filtering  approach  of 
Mitchell  et  al  [3].  However,  the  thickness  of  cloud  cover  certainly 
does  change  over  the  extent  of  an  image,  and,  as  we  stated  in  the  previ¬ 
ous  section,  seems  to  exhibit  rather  deterministic  properties.  Conse¬ 
quently,  an  adaptive  approach  should  be  better  suited  to  restoring  the 
desired  image,  i.e.,  an  approach  where  the  homomorphic  filter  changes 
with  local  cloud  characteristics  which  are  viewed  as  non-stationary  in 
nature. 

Before  proceeding  with  this  adaptive  processing,  in  order  to  make  a 
fair  comparison,  we  shall  first  do  the  best  we  possibly  can  with  the 
stochastic  long-space  approach  of  Mitchell  et  al.  The  improvement  that 
we  achieve  in  their  procedure  relies  on  an  iterative  method  of  estimat¬ 
ing  cloud  statistics  (this  method  in  fact  has  been  suggested  but  was 
never  implemented  [ 3 ] ) . 
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3.1  The  Stochastic  Iterative  Approach 

We  saw  in  Section  2.2  that  a  cloudy  image  can  be  modeled  (after  a  log 
transformation)  as  a  sum  of  noise  and  signal  reflectivity  components 
given  in  (6).  Note  that  this  model  assumes  that  illumination  reaching 
the  ground  hits  objects  uniformly.  We  saw,  however,  that  it  is  not  pre¬ 
cisely  true  due  to  the  different  angular  arrangement  of  objects.  Ne¬ 
vertheless,  this  effect  will  be  assumed  lumped  into  the  function  r(n,m), 
since  as  noted  in  Section  2,  the  resulting  illumination  will  contain  an 
important  high  frequency  component. 

The  objective  of  Section  2.2  was  to  obtain  an  optimal  estimate  of  the 
term  log[L-aLr(n,m) ]  in  the  presence  of  the  "colored  noise"  term 
log[t(n,m)].  We  found  that  we  could  do  this  theoretically  by  applying  a 
Wiener  filter  of  the  form 

H(u, v)  -  ?Mp(u,v)  /  Spp(u,v)  (9) 

-  [Spp(u,v)  -  S^u.v)  -  6(u,v)]  /  Spp(u,v) 

where  the  components  of  (9)  are  defined  in  Section  2.2.  Spp(u,v)  was 
approximated  by  the  squared  magnitude  of  the  Fourier  transform  of  the 
noisy  signal,  and  S^Cu.v)  was  given  by  the  crude  approximation  given  in 
(8). 
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An  alternative  approach  involves  iteratively  updating  the  estimate  of 
the  spectral  density  of  the  clouds.  Since  an  estimate  of  the  reflec¬ 
tance  function  aLr(n,m)  is  required  for  the  cloud  spectral  density 
(see(8>),  we  can  use  an  estimate  of  aLr(n,m)  to  update  the  estimate  of 
the  cloud  spectral  density.  The  Wiener  filter  in  (9)  is  then  updated 
and  the  original  cloudy  image  is  again  filtered.  If  we  do  this  repeat¬ 
edly,  the  resulting  iterative  process  may  converge  to  a  better  estimate 
than  generated  by  the  one  step  process.  In  the  space  domain,  the  ini¬ 
tial  cloud  estimate  is  given  by 

L-s(n,m) 

t.  (n,m)  »  - -  (10. a) 

L-G0 

where  GQ  is  a  constant  (see  Section  2.2).  From  t^n.m),  we  obtain  a 
better  estimate  of  aLr(n,m)  denoted  by  G^n.m),  and  the  iteration  is 
continued  as 

L-s(n,m) 

t.  ..  (n,m)  -  -  (10. b) 

L-Gk(n,m) 

where  Gk(n,m)  is  the  kth  estimate  of  the  function  aLr(n.m).  Figure  5 
outlines  the  procedure.  Figure  6  illustrates  a  modification  of  this  it- 


eration  where  rather  than  filtering  the  original  cloudy  image  on  each 
iteration  we  filter  the  updated  signal  estimate  (i.e.,  filter  G^+1(n,m) 
by  using  Gk(n,m)  to  create  the  new  signal  estimate  G^^Cn.m)).  In  the 
space  domain,  the  iterative  cloud  estimate  is  given  by 

L_Gv+1  <n»m> 

t,  ,  ~  (n , m)  -  -  k“l , . . .  (11. a) 

L-Gk(n,m) 

The  initial  conditions  t^(n,m)  and  t2(n,m)  are  generated  as 
L-s(n,m) 

t  (n.m)  -  (11. b) 

L-G0 

where  GQ  is  a  constant 
and 

L-s(n,m) 

t2(n,m)  -  -  (11. c) 

L-G1(n,m) 

In  theory,  this  procedure  is  not  as  sound  as  that  which  filters  the  ori¬ 
ginal  noisy  signal,  but  as  we  will  see,  it  appears  to  generate  a  better 
restoration.  Since  Gk+^(n,m)  is  a  better  estimate  than  Gk(n,m)  of  the 
original  image,  we  might  use  it  to  obtain  a  better  estimate  of  the  cloud 
transmittance  in  order  to  filter  Gk(n,m)  (i.e.,  swap  Gk+^(n,m)  and 
Gk(n,m)  in  11. a). 
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Fig.  6.  Iterative  procedure  which  operates  on  an  updated  Image  esti¬ 
mate  . 
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3.2  The  Deterministic  Adaptive  Approach 

Consider  now  applying  a  short-space  window  to  the  noisy  signal  P(n,m) 
in  (6).  More  specifically,  let's  assume  the  window  takes  on  a  pyramidal 
shape,  has  a  MxM  square  region  of  support,  and  is  shifted  over  the  data 
at  intervals  of  half  its  width  (i.e.,M/2)  in  both  the  n  and  m  direc¬ 
tions.  If  w(n,m)  denotes  the  pyramidal  window  centered  at  the  origin, 
then  a  particular  windowed  segment  takes  on  the  form: 

Pj  k(n,m)  *  w(n-jM/2,  m-kM/2 )P(n,m)  (12) 

where  the  window  w(n,m)  has  the  following  desirable  property: 

Z  Z  w(n- jM/2 ,  m-kM/2)  -  1  (13) 

.1  k 

Consequently,  the  sectioning  procedure  1b  reversible;  i.e.,  if  no  pro¬ 
cessing  is  applied  to  each  section,  then  the  image  can  be  recovered  ex¬ 
actly. 

Applying  an  adaptive  filter  which  operates  on  each  segment  P  .  (n,m), 

J  »  K 

we  obtain 

Pj  k(n,m)  -  Pj  k(n,m)**hj>k(n,m,a)  (14) 

“  [w(n-jM/2,m-kM/2;P(n,m)]**h  .(n,m,a) 

J  »  K 
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where  the  parameterized  filter  impulse  response  h.  ,  (n,m,0t)  is  a  func- 

J  >  ^ 

tion  of  the  vector  a,  which  is  a  set  of  parameters  dependent  on  the  lo¬ 
cal  cloud  characteristics.  In  particular,  a  is  a  function  of  the  DC 
level  of  the  windowed  signal  which  reflects  the  cloud  density  under  the 
window.  Further,  the  filter  is  high-pass  where  the  shape  and  amplitude- 
depend  on  ou  -The  following  section  will  elaborate  on  this  design.  From 
(6),  the  filtering  process  is  then  given  by 


Pj  k(n,m)  -  [w(n-jM/2,  m-kM/2 ) [ log( t(n,m)+log(L-aLr(n,m) ) ] ] 

**h j  k(n,m,a)  (15) 


We  assume  that  the  window  is  "sufficiently  smooth"  so  that  the  low- 
pass  nature  of  the  noise  term  log[t(n,m)]  and  the  high-pass  nature  of 
the  signal  component  log[L-aLr(n,m) )  are  preserved  after  windowing. 
Consider  now  the  case  where  the  noise  and  the  signal  are  exactly  dis¬ 
joint  In  frequency,  and  the  filter  h,  (n,m,cx)  is  an  ideal  high-pass 
filter  whose  non-zero  energy  band  matches  that  of  the  signal.  Then  we 
can  write 

A 

P  (n,m)  *  w(n- jM/2,  m-kM/2)  log[L-aLr(n,m)]  (16) 

J  > 
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Since  the  noise  and  signal  are  only  approximately  disjoint,  the  noise 

term  will  not  be  entirely  removed,  and  since  the  filter  is  not  an  ideal 

high-pass,  the  signal  will  not  remain  intact.  Alternately,  when 

h  (n,m,a)  adaptively  amplifies  the  high  frequencies  and  attenuates  the 
J  >  k 

low  frequencies,  the  cloud  will  only  be  partially  suppressed  (while  also 
changing  the  low  frequency  component  of  log(L-aLr(n,m) )*)  and  the  high 
frequency  detail  of  the  desired  signal  will  be  enhanced**  (while  boost¬ 
ing  any  cloud  energy  in  this  region).  The  relation  in  (16)  is  then  only 
an  approximation  which  represents  an  estimate  of  the  desired  windowed 
signal  (or  desired  enhanced  signal)  along  with  the  effects  of  any  resi¬ 
dual  cloud  noise.  The  parameters  u  ,  in  h,  (n,m,a),  which  rely  on  the 
local  cloud  density  should  be  chosen  so  that  the  noise  is  suppressed  at 
much  as  possible,  while  the  estimate  of  the  undegraded  image  is  optimal 
in  some  sense. 


*  We  assume,  however,  that  these  components  are  not  perceptually  impor¬ 
tant. 

**  As  shown  in  the  following  section,  under  a  low  contrast  condition 
(i.e.,  with  thick  clouds)  and  with  a  filter  magnitude  greater  than  unity 
in  high-frequency  regions,  the  local  constrast  will  be  enhanced. 


22 


With  these  approximations  In  mind,  our  reconstruction  procedure  which 
we  refer  to  as  overlap-add  can  then,  with  the  use  of  (13),  be  written  as 


P(n,m)  »  £  £  w(n-Mj/2,  m-Mk/2)  log(L-aLr(n,m) ) 

3  k 

*  log[L-aLr(n,m) ]  £  2  w(n-Mj/2,  m-Mk/2)  (17) 

3  k 

=  log(L-aLr(n,m)] 

This  procedure  is  illustrated  in  Figure  7.  Finally,  to  obtain  an  esti¬ 
mate  of  the  signal,  we  exponentiate  and  recover  aLr(n,m)  from 
log[L-aLr(n,m)] . 


W1T(n,m)  OC, 


ln???o-N| 


WLT  -  TRIANGULAR  WINDOW 
WLR  -  RECTANGULAR  WINDOW 


Fig.  7.  Image  restoration  by  deterministic  homomorphic  filtering. 


i 
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3.3  Filter  Design 

From  the  simple  model  of  imaging  above  clouds,  we  have  concluded  that 
it  is  desirable  to  adaptively  attenuate  the  low  frequencies  and  to  am¬ 
plify  the  high-frequencies  of  log  [L-s(n,m)]  given  in  (6).  In  order  to 
proceed  with  the  high-pass  filter  design,  let's  first  look  more  careful¬ 
ly  at  the  process  of  filtering  (6)  adaptively.  To  simplify  the  compli¬ 
cations  due  to  the  logarithmic  operation,  we  shall  approximate  (6)  by 
its  Taylor  series  expansion.  This  is  an  accurate  approximation  under  a 
low  contrast  condition.  In  addition,  we  shall  extend  the  model  in  (6) 
by  allowing  the  downward  attenuation  factor  a  to  take  on  a  space-varying 
characteristic  of  the  form  a(n,m).  This  generalization  accounts  for  the 
possible  changes  in  illumination  reaching  the  ground  due  to  the  varying 
atmospheric  and  other  effects  from  the  clouds  to  the  ground.  Thus,  we 
account  for  a  varying  cloud  thickness  through  t(n,m)  and  other  nonsta¬ 
tionary  effects  (l.e.,  fog,  haze,  etc.)  through  a(n,m). 

The  windowed  segment  in  (15)  then  becomes  : 

Pj  ^(n.m)  -  Wj  k(n,m)  [log[t(n,m>  ] 

+  log(L]+  a(n,m)r (n,m)  ]  (18) 
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If  t(n,m)  and  a(n,m)  are  slowly-varying  over  the  window  w  (n,m),  then 

J » k 

in  the  frequency  domain  our  filter  should  decrease  the  large  negative 
low  frequencies  of  the  Fourier  transform  of  log[t(n,m)]  and  should  in¬ 
crease  the  highs  because  of  the  attenuation  of  r(n,m)  which  holds  per¬ 
ceptually  important  high-frequency  information.  Returning  to  (4),  we 
see  that  we  have  effectively  enhanced  the  local  contrast  in  areas  of 
thick  clouds  by  reducing  the  shift  in  the  local  brightness  level  from 
light  bouncing  off  clouds,  and  by  compensating  for  local  attenuation 
from  passage  of  light  through  clouds. 

Our  filter,  therefore,  should  be  high-pass  with  three  specific  desi¬ 
rable  characteristics.  First,  the  high-frequency  end  of  the  filter 
should  increase  monotonically  as  the  cloud  thickness  increases.  Furth¬ 
ermore,  we  may  wish  to  overcompensate  for  high-frequency  loss  (i.e.,  en¬ 
hance  rather  than  simply  restore  the  image)  by  boosting  the  highs  beyond 
what  would  occur  in  an  undegraded  (cloudless)  image.  Secondly,  we  want 
also  to  increase  the  low-frequency  end  with  increasing  cloud  thickness. 
This  will,  in  effect,  decrease  the  low-frequency  negative  contribution 
due  to  the  log[t(n,m)]  component.  We  note  that  it  may  be  desirable  to 


alter  the  high-pass  filter  so  that  there  exists  a  narrow  pulse  at  DC  to 


accomplish  the  low-frequency  boost*#  We  found,  however,  that  this  fil¬ 
ter  characteristic  was  not  necessary  to  achieve  high-quality  results. 
Finally,  we  may  also  wish  to  broaden  the  bandwidth  of  the  filter  in  the 
high-frequency  range  as  the  cloud  thickens  to  compensate  for  other  forms 
of  high-frequency  loss,  as  for  example  blurring  due  to  the  scattering  of 
light  through  clouds.  A  filter  which  has  the  desirable  shape  and  is 
computationally  efficient  is  a  circularly  symmetric  Gaussian  which  is 
shifted  to  (u,v)-(  tt,tt  ).  A  cross  section  of  this  filter  is  shown  in 
Figure  8  and  the  filter  is  given  by 


H(u,v)  -  A  exp 


u)  =  vu  +  v 


Fig.  8.  Cross  section  of  the  desired  high-pass  filter. 


*  The  resulting  bimodal  filter  then  resembles  that  derived  through  the 
stochastic  approach  of  Mitchell  et  al  [3]. 


The  filter  size  was  chosen  to  be  fixed  at  16x16.  The  filter  parameters 
A,  B,  and  C  are  functions  of  the  value  of  the  Fourier  transform  of  the 
windowed  data  at  DC,  i.e.,  at  (u, v)=*(0 , 0) ,  which  is  an  indication  of  the 
cloud  thickness.  To  determine  the  parameters,  we  first  specify  the  va¬ 
lues  at  (u,v)  =  (0,0)  and  at  (u ,v)  =  (t; ,t  )  for  DC*0  and  for  DC=255.  Each  of 
the  two  pairs  (i.e.,  one  pair  for  the  value  at  (u,v)=(0,0)  for  DC=0  and 
DC”255  and  the  other  pair  than  for  the  value  at  (u,  v)  =  (fT  )  for  DC=0  and 
DC*255  ),  is  then  quadratically  interpolated  according  to  the  formula 

Y  »  (Y2  -  Yr)  D2  /  2552  +  Yx  (20) 

where  D  denotes  the  DC  level  which  takes  on  values  from  0  to  255  and 
where  Y£  and  Y^  denote  two  extreme  values,  one  at  (u,v)  =  (0,0)  and  the 
other  at  (u,v)-(tf , tt ) .  This  yields  the  filter  endpoints  for  DC  levels 
between  0  and  255.  A  quadratic  interpolation  for  computing  the  desired 
intermediate  values  was  used  rather  than  a  linear  interpolation  because 
a  larger  range  of  low  level  of  luminance  needs  almost  no  processing. 
The  filter  parameters  A  and  C  in  (19)  are  then  chosen  so  that  the  filter 
meets  these  computed  endpoints.  The  parameter  B  in  (19)  is  computed  di¬ 
rectly  in  terms  of  the  quadratic  interpolation  formula  (20)  and  two  de¬ 
sired  extreme  values. 
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4.  EXPERIMENTAL  RESULTS 


4.1  Stochastic  Approach  (original  +  iterative) 

This  experiment  was  performed  in  order  to  compare  the  stochastic  ap¬ 
proach  to  the  various  deterministic  procedures.  Figure  9  is  the  result 
of  the  algorithm  of  Mitchell  et  al  where  the  estimate  of  the  typical 
ground  reflection  G(n,m)  was  chosen  to  be  a  constant  equal  to  140.  The 
processed  image  seems  to  contain  primarily  the  high  frequencies  of  the 
unprocessed  image. 

By  using  the  iterative  approach  of  (10),  one  can  improve  the  estimate 
of  the  noise  spectrum.  Figure  10. a  contains  the  result  after  3  itera* 
tions.  Figure  10. b  illustrates  the  result  after  three  iterations  which 
improve  the  estimate  of  the  noise  while  processing  the  image  estimate 
itself  (see  the  iterative  procedure  given  in  (11)).  One  can  see  that 
the  second  form  of  iteration  appears  to  yield  a  perceptually  more  desi¬ 
rable  result.  Finally,  Figure  10. c  depicts  the  result  after  using  *•  it- 
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4.2  Long  Space  Homomorphic  Filtering  (Deterministic) 

In  this  section,  long-space  homomorphic  filtering  was  examined  with 
two  different  filters.  Each  has  a  Gaussian  shape  as  described  In  Sec¬ 
tion  3.2.  The  parameters  of  the  first  filter  are  A“3.65,  B*320,  and 
0-2.15  and  its  shape  Is  given  in  Figure  11.  The  result  is  illustrated 
in  Figure  12. 


Fig.  11.  Cross  section  of  the  filter 
used  in  the  long  space  homomorphic 
filtering 


Fig.  12.  Result  of  the  long  space 
homomorphic  filtering  by  using 
Hj  (u>)  in  Fig.  11 . 


i 
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The  value  at  (u,v)»(0,0)  for  the  two  filters  is  less  than  one  alt¬ 
hough  theoretically  it  should  be  greater  than  one  as  in  the  adaptive  ho¬ 
momorphic  approach.  Since  we  apply  the  filter  on  a  long  space  basis  and 
rescale  the  result  into  the  range  of  0-255,  only  the  relation  of  H(0,0) 
to  H(  tt  ,  tt  )  is  important  as  illustrated  in  Figure  15. 


1 112264-1 1 

RANGE  OF  GRAY  LEVELS  - 

IN  THE  ORIGINAL  IMAGE 
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RANGE  OF  VALUES  IN  THE 
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WITH  H(O.O)  <1 


/ 

/ 
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/ 

/  RANGE  OF  VALUES 
i  AFTER  RESCALING 


255 
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Fig.  15.  An  example  that  illustrates  the  importance  of  the  relation 
between  H(0,0)  and  H f tt  ,tt ) . 
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4.3  Adaptive  Homomorphic  Filtering 


In  this  approach,  a  pyramidal  window  of  size  16x16  was  chosen  and  the 
range  of  parameters  is  given  in  Table  I.  Table  I  requires  some  elabora¬ 
tion,  First,  we  choose  to  allow  the  value  of  the  filter  at  (u,v)=(0,0) 
to  take  on  values  less  than  unity  in  order  to  brighten  extremely  dark 
areas  of  the  specific  image  which  we  are  processing.  Recall  from  Sec¬ 
tion  3.3  that  the  signal  is  log[L-S (n,m) ]  (i.e.,  bright  regions  are 
transformed  into  dark  regions  and  vice  versa).  In  such  a  case  an  in¬ 
crease  in  the  value  of  the  high-pass  filter  at  (0,0)  lowers  the  result¬ 
ing  mean  level  of  the  image  in  the  space  domain).  Secondly,  in  order 
not  to  over- brighten  the  very  dark  areas  of  the  image  we  have  set  the 
value  of  the  filter  at  (u,v)“(0,0)  to  about  0.8  for  all  values  of  DC 
less  than  170,  thus  not  allowing  the  lower  limit  of  0.53  to  be  reached. 
These  parameter  values  were  obtained  after  some  exhaustive  twiddling  to 
ootain  an  "optimal”  enhancement.  We  should  note  that  the  quality  of  the 
resulting  enhanced  image  appears  to  exhibit  a  large  sensitivity  to  small 
perturbations  in  the  chosen  parameters. 

The  local  filter  was  generated  by  interpolating  between  the  two  ex¬ 
treme  filters  (shown  in  Fig.  16)  according  to  the  interpolation  formula 
given  in  (20).  The  result  after  histogram  equalization  is  illustrated 
in  Figure  17. 
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TABLE  1 


ILiNGE  OF  PARAMETERS  FOR  THE  ADAPTIVE  HOMOMORPHIC  FILTER 


DC 

0 

255 

(7!,TT) 

1.1 

1.5 

(0,0) 

0.53 

1.2 

B 

5 

20 

Fig.  16.  Two  extreme  filter  shapes 
which  define  the  limits  of  the 
local  filters  use  1  in  the  adaptive 
homomorphic  filtering  approach. 


Fig.  17.  The  cloudy  image  processed 
by  the  adaptive  homomorphic  filter 
given  in  Fig.  16. 


4.4  Adaptive  High-Pass  Filtering 


In  the  final  experiment,  we  applied  adaptive  high-pass  filtering  to 
the  windowed  image  s(n,m)  in  (A),  rather  than  the  windowed  logarithm  of 
(I.-s(n,m))  in  (6).  The  hig:  frequency  end  ( i  .e. ,  (u,  v)  =  (  r  ,  :  ))  of  the 
high-pass  filter  was  increa:;.  i  monotomically  with  cloud  thickness,  while 
the  low-end  (i.e.,  (u,v)-(n,  ))  was  decreased  monotonical ly .  It  should 
be  pointed  out  that  these  rameter  choices  were  based  on  the  additive 
model  of  (A)  rather  than  tbt  multiplicative  model  of  (5)*.  The  extreme 
filter  values  (again  chosen  by  twiddling)  and  the  resulting  enhanced  ’-m- 
age  after  histogram  equalisation  are  depicted  in  Table  II  and  Figure  18, 
respectively.  From  Figure  18  one  can  see  the  result  is  sharper  than  the 
result  from  adaptive  homomorphic  filtering  but  exhibits  nure  noise  espe¬ 
cially  in  the  darker  regions. 


*  Nevertheless,  applying  a  nonlinearity  (for  dynamic  range  manipulation) 
in  cascade  with  a  high-pass  filter,  from  Schrieber's  studfes  [6],  we 
suspect  that  under  a  low-contrast  condition,  the  resulting  enhanced  im¬ 
age  should  approach  that  obtained  from  adaptive  homomorphic  filtering. 


4. 5  Conclusions 

In  this  section,  we  have  presented  some  informal  preliminary  experi¬ 
mental  results  to  illustrate  the  various  enhancement  procedures  which 
were  proposed  in  the  earlier  sections  of  this  report.  Inspection  of 
Figures  9  through  17  illustrates  that  the  adaptive  homomorphic  procedure 
does  better  in  removing  light  cloud  cover  (i.e.,  in  exposing  objects  un¬ 
der  light  cloud  cover)  than  does  the  stochastic  iterative  or  long-space 
homomorphic  approach. 

Furthermore,  the  adaptive  homomorphic  procedure  appears  to  yield  a 
less  noisy  result  since  a  very  high-frequency  gain  is  not  applied  over 
the  entire  image.  On  the  other  hand,  the  result  from  the  adaptive  homo¬ 
morphic  procedure  contains  less  high-frequency  detail  and  somewhat  more 
low-frequency  "patchiness”  (probably  due  to  the  sensitivity  of  the  local 
mean  level  arising  from  the  nonlinear  logarithmic  operation)  than  that 
from  the  two  non-adaptive  techniuques. 

When  compared  with  the  adaptive  high-pass,  the  adaptive  homomorphic 
yields  a  less  noisy  result,  but  also  with  less  high-frequency  detail, 
particularly  in  darker  regions.  This  is  probably  due  to  the  fact  that 
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in  adaptive  high-pass  filtering,  the  high-frequency  gain  is  applied  di¬ 
rectly  to  the  image,  and  thus  its  effect  is  not  (potentially)  attenuated 
by  the  presence  of  a  nonlinearity.  Finally,  the  adaptive  homomorphic 
procedure  appears  to  do  as  well  as  the  adaptive  high-pass  algorithm  in 
exposing  objects  under  light  cloud  cover. 


5.  SUMMARY  AND  DISCUSSION 

In  this  report,  we  have  demonstrated  the  use  of  adaptive  homomorphic 
filtering  in  exposing  objects  under  light  cloud  cover.  The  adaptive 
filter  Is  high-pass  in  nature  and  Is  parameterized  by  the  local  mean 
level  of  the  degraded  image  which  serves  as  an  indication  of  the  extent 
of  the  local  cloud  cover. 

We  have  shown  through  some  preliminary  experiments  that  this  approach 
compares  favorably  with  the  long-space  stochastic  (iterative)  and  homo¬ 
morphic  approaches.  More  specifically,  the  adaptive  homomorphic  proce¬ 
dure  appears  to  have  a  greater  potential  than  the  long-space  methods  in 
exposing  objects  beneath  light  cloud  cover  while  it  does  equally  well  as 
the  adaptive  high-pass  technique.  Furthermore,  while  the  adaptive  homo- 
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morphie  scheme  does  well  In  suppressing  high-frequency  noise,  it  exhi¬ 
bited  less  high-frequency  detail  than  the  additive  high-pass  filtering 


procedure . 

This  report  has  generated  a  number  of  additional  questions.  Specifi¬ 
cally,  developing  methods  to  suppress  noise  while  exposing  objects  (and 
enhancing  the  detail  of  these  objects)  is  of  paramount  importance.  The 
adaptive  homomorphic  scheme  appears  to  have  a  greater  potential  in  sup¬ 
pressing  noise  than  for  example  adaptive  high-pass  filtering.  One  pos¬ 
sible  reason  for  this  is  that  adaptive  homomorphic  filtering  performs 
local  contrast  enhancement  on  a  nonlinear  function  of  the  image  while 
the  adaptive  high-pass  filtering  can  be  approximated  by  a  differentia¬ 
tion  operation  applied  directly  to  the  Image. 

Another  important  question  involves  the  way  in  which  the  image  is 
segmented.  We  have  in  this  report  divided  the  image  into  overlapping 
square  segments.  An  alternate  method  is  to  process  the  image  adaptively 
on  a  point-by-point  basis  (similar  to  the  procedure  of  Peli  and  Lim 
[4]).  It  is  important  to  understand  the  differences  in  these  approach¬ 
es,  with  respect  to  their  computational  and  noise  suppression  character¬ 
istics. 
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